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ERROR ANALYSIS OF HOUSEHOLDER TRANSFORMATIONS AS APPLIED 
TO THE STANDARD AND GENERALIZED EIGENVALUE PROBLEMS 


By Robert C. Ward 
Langley Research Center 

SUMMARY 

Backward error analyses of the application of Householder transformations to both 
the standard and the generalized eigenvalue problems are presented. The analysis for 
the standard eigenvalue problem determines the error from the application of an exact 
similarity transformation, and the analysis for the generalized eigenvalue problem 
determines the error from the application of an exact equivalence transformation. 

Bounds for the norms of the resulting perturbation matrices are presented and compared 
with existing bounds when known. 


INTRODUCTION 

Examination of the eigenvalue algorithms which are recommended by Wilkinson 
and Reinsch (ref. I) for solving the various classes of eigenvalue problems reveals that 
Householder transformations, sometimes referred to as elementary unitary Hermitian 
transformations (see ref. 2), are used extensively in these algorithms. A Householder 
transformation can be represented by the expression I - A wT where I is the iden- 
tity matrix, v is a vector, and c is a scalar equal to (l/2)v T v. These transforma- 
tions are normally used to transform a vector x (usually a portion of a column of a 
matrix) into a vector y which has only one nonzero component (usually the first 
component). 

Wilkinson (refs. 2 and 3) and Ortega (ref. 4) have published detailed error analy- 
ses of the application of Householder transformations in the standard eigenvalue prob- 
lem based on the same general approach of determining the error from the application 
of an exact unitary similarity transformation. These error analyses are important and 
have earned their reputation in numerical linear algebra. Since eigenvalues are pre- 
served by similarity transformations whether these transformations are unitary or not, 
another realistic approach is one of determining the error from the application of an 
exact similarity transformation. 



In this paper, this latter approach is examined and a backward error analysis is 
presented first for Householder transformations in the standard eigenvalue problem. 
Then, a backward error analysis is presented for these transformations using the same 
approach but applied to the generalized eigenvalue problem. 

SYMBOLS 

When one symbol is related to another symbol, this relationship is identified 
clearly by the context in which the symbol appears, and although some symbols have 
multiple definitions, the context makes the meaning of the symbol unambiguous. 

A,B,D n x n real matrices 

n x n perturbation matrices 
I identity matrix 

P,Q,Z n x n Householder transformation matrices 
S a scalar used to compute a Householder transformation 

c scalar in a Householder transformation 

i,j,k,r non-negative integers 

n order (size) of the matrices 

p, q scalars used in lemma 1 and lemma 2 

q, u,v n x 1 vectors 

i number of binary digits used to represent the mantissa of a floating-point 

number in a computer 

6q,6u,w n x l error vectors 

e ,V scalars representing the relative error in a computation 


e»f,g,h,'i 

X,Y,Z j 
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an eigenvalue 


\ 


Superscripts: 

-1 inverse 

i exponentiation to ith power 

T transpose 


Subscripts: 

i ith item in a sequence or ith component of a vector 

ij element in (i,j) position of a matrix 


Special notation: 

0( ) order of magnitude 

| [ absolute value (if V is a vector, |V[ is the vector with components |V t 

If A is a matrix, |A| is the matrix with elements j ^ij j* j 

|| | | unspecified norm 

j| || 2 two norm (if V is a vector, ilvllg is the value ^ |Vi| 2 j 1//2 . If A is 

a matrix, ||A|| 2 is the value (max \^V 2 where \ are the eigen- 
values of a?a) 


II 


ll E 


Euclidean norm of a matrix 



is the value 



replacement symbol (computer equal sign) 

A tilde (~) over a matrix indicates an n x n matrix related to that matrix. 
A prime (*) over a matrix indicates an n x n matrix related to that matrix. 



A horizontal bar (”) over a symbol indicates the computed value, whereas an exact 
number is indicated by a symbol without the bar. 

Other mathematical notation has its usual meaning. 

PRELIMINARIES 


A basic assumption which is used in the error analysis is the following inequality: 

n2-t < 0.006 (1) 

where n is the size of the matrix and t is the number of binary digits used in rep- 
resenting the mantissa of a floating-point number in the computer. This assumption 
compares with that of n2 _ t < 0.00032 used by Wilkinson in reference 3 and 
n2 -t < 0.008 used by Wilkinson in reference 2. Under the assumption given by equa- 
tion (1), n would be restricted to be less than 1.68 x lO 1 ^ 0 n the Control Data series 
6000 computer systems which have 48 binary digits in the mantissa of a floating-point 
number. In most scientific computers, the restriction imposed by limited computer 
memory therefore automatically assures that equation (1) will be satisfied. 

In addition, the following inequality is assumed to be valid: 

2-t<2~ n ( 2 ) 

or t s ii. This assumption compares with that of 2 _ * < used by Wilkinson 

(ref. 3) and 2“* = 10 used by Ortega (ref. 4). 

The following two lemmas are used extensively in the analyses. 


Lemma 1: If 0 < p < 16, then 

(l + p2 _t ) n < 1 + (1.06)np2- t 

Proof: Let q = p2 _ *. 

By the hypothesis 0 < p, then 0 < q and 


1 + q < e^ = 


1 


+ q + 


q£ 

2 ! 



+ . . . 


Therefore, 

(1 + q) n < e«l 

Thus, 


(1 + q) n < 1 


+ nq 



l »q J 


( 3 ) 
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If the exponential expansion is used, then 


e n q=1 + nq + (^! + M! + 


and 

e nc i - 1 _ i , , (°q ) 2 , ( n q ) 3 , 

nq 21 31 4! 

Comparing these two expansions term by term and using the hypothesis that 0 < p or 
0 < q, yields 


or 


g nq -1 < l(e n q + 1 ) 
nq 2 


Thus, from equation (3) 

(1 + q) n < 1 + nql 


i(e n( l + l) 


( 4 ) 


(5) 


When the assumption given in equation (1) and the hypothesis p < 16 are used, then 
nq < 0.1 

and 

I(e nq + l) < 1.0526 (6) 

Combining equations (5) and (6) and the definition of q produces the desired result. 
Lemma 2: If 0 < p < 16, then 

(l -p2' t ) n > 1 - (1.06)np2- t 
Proof: Let q = p2 _t . 

By the assumption given in equation (2) and the hypothesis 0 < p < 16, then 
0 < q < 1 and the following expansion is valid: 

(1 - 1 )° = 1 - nq 4 . q2 - n(n ~ » (rl ~ 2 > q3 + . . . „ q n-l ± q „ 


3! 


> 1 . nq . . Mi . . . . . Mid . Mi > 2 . e "q 

21 3! (n - 1): n! 
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Thus, 


(1 - q) n > 1 - nq 


e 


m . l 
nq _ 


Since equation (4) is based only on the hypothesis of 0 < p or 0 < q, then this equation 
is also valid here. Hence, 


(1 - q) n > 1 - nq 


i(e"<I * 


1 


( 7 ) 


When the assumption given in equation (1) and the hypothesis p < 16 are used, then 
nq < 0.1 

and 

I(e n( l + l) < 1.0526 (8) 

a 

Combining equations (7) and (8) yields the desired result. 


ERROR ANALYSIS OF HOUSEHOLDER TRANSFORMATIONS 
IN THE STANDARD EIGENVALUE PROBLEM 


The error analysis for Householder transformations in the standard eigenvalue 
problem Ax = \x is subdivided into the analysis of three errors. First, the error 
from computing a similarity transformation under exact matrix multiplication is pre- 
sented. Then, the error in the matrix multiplication is analyzed. Finally, these results 
are combined and an error analysis of a sequence of Householder similarity transforma- 
tions is presented. 


Error Analysis of Similarity Transformation 

Let a computed Householder transformation be given in terms of a scalar c 
and a vector v by 

Pi = 1 . i vv T (9) 

(Details concerning the derivation of Householder transformations are found in ref, 2.) 
Because Householder transformations are theoretically unitary and hermitian, similarity 
transformations are performed by premultiplying and postmultiplying a matrix by Pj 
as defined by equation (9). 

Consider the problem of applying the Householder transformation Pj to A^.j 
as a similarity transformation and denoting the resulting matrix by A^ where the 


6 



matrix multiplications are computed exactly; that is, 


A| - Pi^i-iPi 

This expression will be a similarity transformation if Pj = Pf Let 

Pi = Pf 1 + Ei (10) 

Then 

A t = (Pf 1 + E^Ai.jPi = Pf 1 (A i .i + PiEiAi.^Pi 

If a bound for |PiEjAi_iJ| could be found, then this would be a bound on the perturba- 
tion added to Ai_i in order to make A^ _ ^ and A^ similar matrices. 

Theorem 1: Let e be the relative error in computing c given the vector v; 
that is, 

c = c(l + e) 

where 

c = -i v T v 


Then 

||P i |l 2 ^l + 2|ei + 0(e 2 ) 

Proof: Since Pj is hermitian, j|Pi|| 2 = maxj\jj where are the eigenvalues 

of Pi. The Householder transformation Pi has an eigenvalue equal to 1 with multi- 
plicity of n - 1 and the remaining eigenvalue is 1-4 v^v. Thus, 



max 



Using the hypothesis c = c(l + e) with c = -i v^v and the fact that the transformation 
would not be applied if v T v = 0 yields 



max 




= max 



5 1 + 2 Id + o(e 2 ) 


Theorem 2: Let e be as given in theorem 1. Then 
||Pf 1 || 2 = 1 + 2 |e| + 0(e2) 
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Proof: Since Pf 1 is hermitian, Pf 1 ) 2 = maxjX-jJ where \j are the eigen- 
values of Pj-1. Also the eigenvalues of Pr-1 are the reciprocals of the eigenvalues 


of Pj. Thus, 


Ml 2 = max ll, (l - 4 vTv)' 1 


Using the hypothesis c - c(l + e) with c = A v T v and the fact that the transformation 
would not be applied if v^v = 0 produces 


1 + e> 


1 < 1 

’!! - 6 / 


+ 2[e| + o(e 2 ) 


Note that a formal expression for Pj" 1 can be given; that is, 


P f A = I + 


1 --T 

± ir\7 -I- 


c - V T V 


Theorem 3: Let e be as given in theorem 1. Then 
j|E i iL5 4| C | + 0(|e|3) 


Proof: From equations (9), (10), and (11), 

iWi E = llpr 1 - Pill E 


I + W T -1+4 VV T 

C - V^V c j; 


+ 4 IIvvT|I e 

c - v x v c 


~ + 4 v T v 

c - V i V c 


Because of the definitions of e and c, the norm may be written as 


1 E 1 -T- 

r v 


i v T v(l + e) - v T v A v T v(l + e) 


+ v^v = ? h — — 


-1 + 6 1 + € I + e 2 


£ 4161 + OUel 


Note that Parlett (ref. 5) has presented a similar result to theorem 3 in a slightly differ- 
ent context. 
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Thus, as theorem 3 shows, a bound for the norm of the perturbation resulting 
from one Householder similarity transformation with exact matrix multiplication can 
be exhibited. Since Pj is a hermitian and normal matrix, then ||PjX| 5 ||Pj|L||x[|g 
for any matrix X, and 

' = [l + 2 |e| + o(e 2 )][4[e| + -0(|el 3 )j ||Ai_i|| E 
?4|e| ||Ai_i|| E + o(|e|2) 

For small e, an approximate bound can be given; that is, 

jPiEiA U1 j! E r 4.6: |jA ul || E (12) 

All that is needed now is a bound on e; that is, a bound on the relative error in 
computing c given the vector v. Since the vector v differs from a vector x 
already in the computer by only the first component, the following algorithm, suggested 
by Parlett (ref. 5), for computing c should be used to obtain a small relative error e: 


(1) 

S - 

- x 2 2 + x 3 ^ + • • 

• +x n 2 


If 

S = 0, skip this 

transformation. 

(2) 

|k| 

- sqrt^Xj 2 + s) 


(3) 

k - 

- - jkj sign(xj) 


(4) 

*1 

- - k 


(5) 

c - 

■ \h 2 * s) 



It is assumed that the computations are made with a 2t -digit accumulator (Wilkinson in 
refs. 2 and 6 uses the notation fl 2 for such computations) and that S computed in 
step (1) retains its 2t digits for use in step (5). 

From the information in reference 6 and the fact that the multiplication of two 
t -digit numbers is exact if a 2t -digit answer is retained, the computed value of S is 
given by the following equation: 

S = x 2 2 (l + e 2 ) + x 3 2 (l + e z ) + . . . + x n 2 (l + e n ) 


where 


(i Jj-arhui,!! 

l + 42-2 ‘) n ' 2 

\ 2 / 2 \ 

2 / 
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and 


(‘-I 2 ' 21 ) 


-2t\ n+1 " r 


z 1 + 


s(l + |2- 2t ) 


n+l-r 


(r=3, 4, 


n) 


Since a bound is desired for the relative error in computing c given the vector v, 
the error in computing vj does not enter in this analysis. Thus, the analysis proceeds 
to step (5) of the algorithm. 

From the information in reference 6 and the foregoing assumptions, 


where 


and 


Vi 2 + S = v x 2 (1 + ej) + x 2 2 (l + e 2 ) + x 3 2 (l + e 3 ) + . . .4 x n 2 (l 4 e n ) 


(l - | 2 _2t j S 1 4 £ (l 4 j 2- 2t ) 

( l -| 2 - 2 ‘) n - 1 S 1 +£ 2 S ( l + | 2 - 2‘) n - 1 


(» -f 2 ' 2 *) 


n42-r 


5 1 4 6 


r S (l + I 2- 21 ) 


n42-r 


(r = 3, 4, 


n) 


A bound on j€j| can be obtained by inspection, and bounds on |e 2 | and |e r j can be 
obtained by applying lemmas 1 and 2 with p = | 2 _t . The bounds are as follows: 



|€ 2 | s|(1.06)(n - l)2' 2t 

l e rl = l[(l-06)( n 4 2 - r)2 -2 ^ (r 5 3, 4, , . n) 

Therefore, 

Vj 2 4 S = (vj 2 4 X 2 2 + x 3 2 4 . . .4 X n 2 )(l 4 e 0 ) 

where 

|e 0 | < |(1.06)n2~ 2t 

Since dividing by 2 in step (5) introduces no error on a binary computer, an expression 
for c can be obtained; that is, 
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c - ^l 2 + x 2 2 + x 3 2 + • • • + x n 2 )(l + e 0 )(l + v) 

where 

|e 0 j < |(1.06)n2- 2t 

and 

; 7 ]| < 2 _t 

The error t] is due to rounding to single precision. 

When the assumption given in equation (1) is used, the error reflected in Cq is 
overshadowed by the error rj. The exact value c, given the vector v, would be equal 
to + xg2 + + . . . + x n 2), and then the relative error in computing c is 

bounded by 2"^; that is, 

Id £ 2-* (13) 

Thus, equation (12) representing a bound for the norm of the perturbation resulting from 
one Householder similarity transformation with exact multiplication can now be given by 

||PlE i A i _ 1 || E £(4)2-t||A i _ 1 || E 

Note that the error bound given for ||p - P”*||g of (4)2 may be compared 
with the error bound given by Wilkinson (ref. 2) for IIP - P'lL of (9.01)2”* : where 
P’ = I - i vv . Then, in a sense, PAP is roughly twice as close to being an exact 
similarity transformation as it is to being an exact unitary similarity transformation. 

Error Analysis of the Matrix Multiplication PiA^iPj 

In the previous section, the errors made in the matrix multiplications PiAi_iPi 
were ignored. In this section, these errors are analyzed. 

Let be the error in the computation of P[Ai_iPi; that is, 

A l = PiAi.jPi + F i (14) 

The error in the premultiplication by Pj is now considered. Wilkinson (ref. 2) has 
presented an error analysis of this premultiplication under a slightly different assump- 
tion for n2 - t as explained previously. This analysis is very similar to that of 
Wilkinson. 

Let P = Pi = I - i vvT, A = Ai_i, and B = PA. (Note that c and v do not 
have bars over them. Since this section is concerned with the matrix multiplications 
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after is formed, the analysis assumes that c and v are exact and not computed 
values with errors.) Then, the premultiplication algorithm becomes 

(1) u T - ^ v t A 

(2) B — A - vu* 1 

It is assumed, as in the previous section, that the computations are made with a 2t -digit 
accumulator. 

When the information in reference 6 and the fact that the multiplication of two 
t-digit numbers is exact if a 2t-digit answer is retained are used, the ith component of 
the computed value of uT in step (1) is given by the following equation: 

G i T = 3pi a li(l + e l) + v 2 a 2i( 1 + e 2) + ■ • • + v n a ni (l + (l + *o) 

where 

(l - 1 2- 2t ) n ' 1 ! 1 * f[ 5 i'l + | 2- 21 )"" 1 

(l - f 2 - 2 t)" +1 - r Sl + ir s(l + | 2 - 2 t)" +1 - r (r - 2, 3, . . n) 

and 

(l - 2-*) * 1 + e Q ^ (l + 2-t) 

Application of lemmas 1 and 2 with p = 2“t yields a bound for Je-jj and |e r | that 

can be expressed by -|(1.06)(n - l)2" 2t . A bound for je 0 | can be obtained by inspec- 
tion to be 2 _t . Therefore, 

uj T = ui T (l + e 0 ) + wi 

where 

HI = |( 1 - 06 )( n - l)2- 2t (l + 2- t )l[[v 1 |ja li j + jv 2 j|a 2i | + . . . + |v n J |a ni j] 

Hence, 



and 

llu - uj| 2 s 2- t i|u|| 2 +|(1.06)(n - l)2- 2t (l + 2-t)i||v|i 2 !|a|| E 
With the assumption given by equation (1), the last expression becomes 

IM Z = llu - u(| 2 s 2-t||u|| 2 + (0.01)2-* i||v|l 2 IKIIe ( 15 ) 
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The (i,j) component- of B in step (2) is given by the following equation: 
by = [ay (1 + Vl ) - YiUj (! + >? 2 )j (1 + %) . 

where 

(l - | 2" 2t ) r ^ 1 + 77 r ^ (l + | (r = 1, 2) 

(l - 2-*) s i + Vo < (l + 2-*) 

A bound for and J77 2 | is then | 2" 2t and a bound for |?7 q| is 2 _t . Therefore, 

b ij = b Lj t 1 + ^o) + [ a U ?7 i - v i u J^2 " V i 5u J t 1 + ^2) t 1 + ^0) 

Thus, 

b ij - b ij| = 2_t | b ij| + [| 2 " 2t | a ij| + | 2 ' 2t | v ii) u j| + l v i| H (* + | 2' 2t ) (l + 2-t) 

|§ - B| * 2- t jB| + jj 2- 2t |Al + | 2- 2t |vlM + |v||6u T | (l + | 2' 2t j (l + 2"t) 

and 

I!b - b|I e i 2-tj|B|l E + [j 2-2t|| A || E + | 2-2t||v|l 2 l|u|[ 2 + l!v|j 2 |M 2 (l +| 2- 2 *) (l + 2-t) 

This bound would be more convenient and usable if it could be expressed in terms 
of 1 1 A|| e* From the definition of B, theorem 1, and equation (13), a bound for ||b|| e 
in terms of || Aj| E can be obtained by the following sequence: 

IIb|| E = !|pa|| e s||p|| 2 ||A|| e S [l + (2)2-t]||A|| E (16) 

In the proof for theorem 1, Iv^v is bounded by 2 + 2e or 2(1 + 2~t). Therefore, 
the following sequence is valid: 

rivH 2 ||ull 2 ^'I||vil 2 2 !|Ai| E £ i vT v ||a|| e ^ 2(1 + 2-t) f|A|| E 

When equation (15) and the foregoing sequence are used, the following sequence is obtained 
to bound I Ivllg I |Suj) 2 • , . , . 

Ilvll 2 ||6u|i 2 5 I|vll 2 j2 _t |fui!2 + (0.01)2 _t I ||vll 2 |fA(|.jgj 

? (2^) (2) (l + 2-t)|jA|| E + (0.01) (2 -t-)(2) (l + 2-t)||A|| E 5 (2.03)2-t||A|| E 
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A bound on IIb - Bllg can now be obtained in terms of ||A[I e 5 that IS > 



5 (3.04)2-t||A|| E (17) 

This bound compares with Wilkinson's bound (ref. 2) of (3.35)2-*) |A| under a slightly 
less stringent assumption for n2 - ^. 

The error analysis for the postmultipli cation by Pi is very similar to the pre- 
multiplication. The steps in the postmultipli cation algorithm are 

(1) q - | Bv 

(2) Aj — B - qyT 

If 2t-digit accumulators and lemmas 1 and 2 are used, the ith component of the computed 
value of q is given by 

Si = £Pil*l(l + e l) + S i2 v 2( 1 + e 2 ) + ’ • * + 6 in v n(l + % 

where 

|e r | <|(1.06)(n - l)2-2t (r = 1, 2, . . n) 

and 

N - 2 "‘ 

Hence, 

Si = ^i? + e o) + 

where 

M * |»-06)(n - l)2-2t(l 4 2-t)i(j6 ll ||v 1 | 4 |b 12 ||v 2 | + . . . + |6 ln ||v„] 

Thus, 

1MI 2 5 119 - 4|b S 2 -‘IW| 2 + to- 01 ) 2 ' 1 i||v|| 2 l|B|| E ( 18 ) 

If D represents the exact product BPj and D represents the computed value 
BPi, the (i,j) component of D is given by the following equation: 
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where 


dij = [bij (l + ijj) - qi v j (! + V 2 )\ (! + 


D 


W < 2_t 
N 5 1 2 ‘ 2 * 

M s I 2-21 


Therefore, 

|d,j - dij| S 2-t|dij| + f| 2-2t|bij| + | 2-2t| qi ||vj| + |6 qi ||vj| (l ♦ | 2-2‘j 


(l + 2-t) 


and 


||6 - D|! E 5 2~ t ||D| | E + (| 2- 2t |}B|| E + | 2- 2t ||q|| 2 l|v|f 2 + ||Sq|[ 2 l|v|| 2 (l + | 2' 2t ) (l + 2"*) 


As in the premultiplication, this bound would be more convenient and usable if it 
could be expressed in terms of ||b(| e which could then be expressed in terms of ||A|I E . 
From the definition of D, theorem 1, and equation (13), a bound for ||D|| E in terms of 
||bI1 e can be obtained as follows: 

IIDIIe - IIbpIIe § I|b|! e ||P|I 2 ^ [l + (2)2- t ]ltB|| E 


Again, using the bound for I v^v yields the following sequence: 

!|q|t 2 ll v !l2 = ^livl I 2 2 I|bI1 e £ i v T vllB|| E ? 2(l + 2- t )llBit E 
Application of equation (18) and the foregoing sequence leads to a bound for ||5q|| 2 )lvll 2 : 


||6q|| 2 l| V tl 2 ?ll V il 2 


2 _t i |q| U + (0.01)2 -t i| |v| l 2 l IB| |j 


s (2 -t) (2) (l + 2-*) ||b1I e + (0.01)(2- t )(2)(l + 2- t )llBll E 

? (2.03)2' t l!Bl! E 
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Thus, 

IlD - Dllj; 5 

+ (l +| 2-2tj (2.03)2 -^ Hb|I e 

= (3 . 04 ) 2 -t HbI I E (19) 

An expression for the error in the combined premultiplication and postmultiplica- 
tion in terms of ||A|I e can now be obtained from the definitions of B, D, and D 
(eqs. (13), (16), (17), and (19)) and theorem 1. Application of these definitions results 
in the following sequence: 

l| A i - P iVlPi|| E - ||Ai - BPj| E = ||At - BPi + (B - B)Pi|| E 

= U A i - BPi|j E + ||(B - B) Pi |( E £ |jD - D|j E + ||B - Bll E f| Pi [| 2 
* (3.04)2- t llBtl E + (3.04)2- t |!A i _ 1 || E [l + (2)2*3 
£ (3.04)2't||Bll E + (3.05)2-t||A i _ 1 || E 

£ (3.04)2- t jTlB(l E +IIB - Bllgj + (3.05)2- t ||A i _ 1 || E 
£ (3.04)2' t |[l + ( 2 ) 2-9 + (3.04)2“j \\A iA \\ E + (3.05)2- t ||A i _ 1 || E 


^2-t g + (2)2-3 + (l + 2-t)[| 2-2t + | 2-2t(2)(l + 2-*) 


S (3.05)2-t||A l _ 1 |[ E + (3.0S)2-t|fA|. 1 |J E 

Thus, the bound on the error matrix Fj in equation (14) can be given by 

||F 1 || E S(6.1)2-‘||A l . 1 || E (20) 

Error Analysis of a Sequence of Householder Similarity Transformations 

There are several algorithms which require the use of a sequence of Householder 
similarity transformations; for example, the reduction of a matrix to Hessenberg form 
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requires a sequence of n - 2 such transformations. This section combines the results 
presented in the previous two sections to obtain a bound for the norm of the perturbation 
matrix which results from computation of such a sequence. 

Apply a sequence of k Householder transformations to a matrix A. When 
Aq = A, the sequence is denoted by 

Ai - PtAi.iPi (i = 1, 2, . . k) 

For each transformation, 

Ai = (Pi * + Ei'jAi.jPi + Fj = Pj *Ai _iP[ + EiAi_jPi + Fi (21) 

Then, 

A k = p k lp kli • • • + Ej + Fj + E 2 + F 2 + . . . + E k + Fk|PiP 2 . . . P k 

where 

Ei = P i E j^Aq 

Ei = P!P 2 . . . P i Ei Aj _ iP Pf _^2 . . . Pf 1 (i - 2, 3, . . k) 

and 

Fi = PiP 2 . . . PiFiP^Pf. 1 ! . . . Pj 1 (i = 1, 2, . . ., k) 

If the bounds for Pj, Pj" 1 , Ej, and Fi are taken from theorems 1, 2, and 3 and from 
equations (13) and (20), then 

||Ei|| E § [l + (2)2i 2l_1 [(4)2-5||A i .i|| E 


and 


INI e S [l + (2)2-‘] 2, [(6.1)2 -*]||a 1 . 1 || e 

From equation (21) and the bounds for P^ E^, and the following sequence 

results; 

[|Aij| E S Q + (2)2-t) 2 ||A i . 1 || E + [l + (2)2-3 Q4)2-t]||A l . 1 || E + (6.1)2-t||A i . 1 || E 

* j[l + (4)2”* + (4)2‘ 2t ] + |(4)2“ t + (8)2' 2 5 + [(6.1)2^ ||Ai_iH E s fl + (14.11)2-^| |A t _i| | B 
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Therefore, 

INI e S [l + (2)2-‘] 21 - 1 1(4)2-‘] [l - (14.11)2-‘] 1 - 1 ||a 0 || e 

and 

INI e S [l + (2)2' l ] 2i [l + (14.11)2-g l - 1 ||A 0 || E 

Applying lemma 1 yields 

||Ei|| E < [l + (1.06)(2)(2i-l)2- t ] 1(4)2-*] [l + (1.06)(14. ll)(l -1)2"*] 1 1 A 0 | | E 
1(4)2-* + (76.77i)2- 2 * + (253.6i 2 )2' 3 *] ||a 0 | | 

[l + (1.06)(2)(2i)2-*] [(6.1)2-*] [l + (1.06) (14 . 1 1) (i - 1)2 -*]||a 0 |L 

E 

< [(6.1)2-* + (117.1i)2- 2 * + (386.7i 2 )2- 3 *][|A 0 || E 

Note that the i and i* terms have not been deleted since no restriction was placed on 
the number of transformations in the sequence and i may be large enough to prevent 
those terms from being ignored. Hence, 

A k = P" 1 (A + Z)P 

where 

k 

IIZll E = ^ E 10 - 1 ) 2 ' 1 + (193-87)i2- 2 * + (640.3)i 2 2- 3 *] |IAII E 

i=l 

* ((10. 1)k2 "* + (97.0)k 2 2 -2 * + (213.5)k 3 2- 3 *] ||A|I E 

This bound can be illustrated by the reduction of a matrix to Hessenberg form* 
This reduction requires n - 2 Householder transformations* Thus, the bound on the 
perturbation matrix Z which yields the exact similarity transformation is given by 

IIZ|| E £ [(10.1)n2-* + (97.0)n 2 2 -2 * + (213.5)n 3 2- 3 *]||A|| E 

Using the assumption given by equation (1) or n2 _ * < 0.006 gives 
IIZ|I E £ (10.6)n2-*)|A|l E 


< 

and 
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This bound may be compared with a bound of (24.72)n2 _ t||A|| g given by Wilkinson 
(ref. 2). Thus, in a sense, PkPk-1 • • ■ PiAP]P 2 . . . Pk is roughly twice as close to 
being an exact similarity transformation as it is to being an exact unitary similarity 
transformation. 

ERROR ANALYSIS OF HOUSEHOLDER TRANSFORMATIONS 
IN THE GENERALIZED EIGENVALUE PROBLEM 


Because of the development of the QZ algorithm (ref. 7) and the combination shift 
QZ algorithm (ref. 8), Householder transformations are being used extensively in solv- 
ing the generalized eigenvalue problem Ax = \Bx. In this case, the concern is not to 
use a similarity transformation but an equivalence transformation. Thus, there are no 
perturbations on the original A or B matrix due to the premultiplication or postmul- 
tiplication by the Householder transformation Pj instead of Prl; that is, the E^ and 
Ej matrices in the preceding sections are the zero matrix. However, the error in a 
premultiplication given by equation (17) and the error in a postmultiplication given by 
equation (19) are still valid in the generalized eigenvalue problem. 


Apply to A and B a sequence of k premultiplying Householder transformations, 
denoted by Q*, and k postmultiplying Householder transformations, denoted by Zj. If 
Aq = A and Bq = B, the sequence is denoted by 



(i = 1, 2, . . ., k) 


Both the combination shift QZ algorithm and the QZ algorithm proceed in this fash- 
ion. For each Qi and Zj transformation, 

Ai = (QiAi.! + Gj)z t + H| = QiAi.iZj + G^ + H[ (22) 

and 

§i = (QiB ul + Gi') Zi + H^ = QiBi.iZi + Gt’Zi + (23) 


where Gi and Gi' are premultiplication error matrices and Hi and Hi' are post- 
multiplication error matrices. Then 

Ak= QkQk-1 ■ - • Ql|^o + ^1 + **1 + ^2 + + • ♦ • + Gk + HkjzjZ2 . . . Z k 

and 

Bk = QkQk-1 • • • Ql [po + ^l' + ^1* + ^2* + ^2* + • • • + Gk' + Hk'JzjZg . . . 
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where 


Gl = 


Gi = Qj/Q 2* • • • Qj[ ^GiZ^ 3i^i^2 • • • (i — 2, 3, . . k) 

Hi = Q^Qa 1 . . . QiVi^Zii . . • z 'i 1 (i = 1, 2, . . k) 

Gl’ = Qi 1 G 1 ' 


Gi’ = Qi'Qj 1 • 

• • Qf^i’ZfiiZfi • • • z ; 1 

H 

i° 

CO 

. ., k) 

Hi' = . 

■ • Qi^Hi'Zf'Ziij . . . zji 

II 

• k) 


By taking the bounds for Qj" 1 , Zj 1 , G[, and Hi from theorem 2, and using equa- 
tions (13), (17), and (19), one obtains 

INI e S [l + <2)2- t ] 2l ' 1 [(3.04)2-»]||A t . 1 || E 

INI e s 0 + (2)2-‘] 2 ‘|(3.04)2- t ]||Q l A 1 . 1 + 0l || E 

S [l + ( 2 ) 2 -‘] 2 '|( 3 . 04 ) 2 - t ]|j|Q i |y|A 1 . 1 || E+ ||G 1 || E j 


5 [l + (2)2- t ] 2 ‘[(3.04)2- t ] |l + (2)2-'] + (3.04)2 -*J||A,_i 


E 


'E 


§ [l + (2)2-t] 21 [(3.05)2-*] | IViH 
From equation (22) and the bounds for Qj, Z- L , Gj, and Hi, 

[l + (2)2- t ] 2 ||A i _i|| E + [l + (2)2 _t J [(3.04)2- t ]|jA i _i|| E + ((3.05)2-*] 


A i|lE = 


^i-1 


^ [l + (10.1)2~ t ]||A i _ 1 | 


E 
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Therefore, 


ll°i|l B * [i + W2-9 21 ' 1 is.o4)a-9[i + (io.i)a-g iJ ||A 0 || B 

||Hj! E £ [l + (2)2-‘] 2l [(3.05)2-‘] [l + (J0.1)2- l ] 1 " 1 ||A 0 || E 
Applying lemma 1 yields 

INI e < [l + (1.06)(2)(2i - l)2-t] [(3.04)2-*] [l + (1.06)(10.1)(i - 1) 2 |Ao| [ E 

< |(3.04)2 - * + (45.451)2-2* + ( 138 . 10i2) 2 -3t| j | a 0 | ( E 

||Hi|| E < [l + (1.06)(2)(2i)2- t ] [(3.05)2-*] [l + (1.06)(10.1)(i - 1)2 '*] | |Ao|| b 

< [(3.05)2"* + (45. 59i) 2 _2t + (l38.51i 2 ) 2“ 3t ]| | a q | | 

If Q k Q k _i . . . Qj is denoted by Q and . . . Z k is denoted by Z, then 

A k = Q( A + X)Z 

where 


k^ 

l|X|| E ? y [(6.09)2-* + (91.04)i2-2t + (276.61)i22-3*] ||A|( E 
i=l 


£ (6.1)k2-* + (45.6)k 2 2- 2 * + (92.3)k 3 2- 3 *] IIAII 


E 


Similarly, 


where 


B k = Q(B + Y)Z 


||YH e ^ |(6.1)k2“ t + (45.6)k 2 2 _2t + (92.3)k 3 2- 3 *] I |BI( 


E 


SUMMARY OF RESULTS 


Householder transformations have been analyzed with the goal of obtaining new 
bounds on the perturbation matrix in a backward error analysis for both the standard and 
generalized eigenvalue problems. The important results of this study are as follows: 
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1. A bound for the norm of the perturbation matrix was obtained for one Householder 
similarity transformation with exact matrix multiplication. That is, if P is a computed 
Householder transformation and A = PAP exactly, then A = P-1(A + E)P where 

IIEII E 5 (4)2- t llAII E 

2. A bound for the norm of the perturbation matrix was obtained for a sequence of 
k Householder similarity transformations. That is, if Pi,P2> • • • > p k are computed 
Householder transformations, A is the computed P^P^.j . . . PjAP^^ • • • ^k’ anc * 

P is defined to be the exact product P 1 P 2 • • • Pk> ^ en ^ = + F)P where 

IIFII e i |(10.1)k2 -t + (97.0)k 2 2- 2t + (213.5)k 3 2- 3t ]llA|l E 

3. A bound for the norm of the perturbation matrices was obtained for a sequence 
of k premultiplying and postmultiplying Householder transformations with regard to the 
generalized eigenvalue problem. That is, if Q-pC^, • . . , Q k ,ZpZ 2 , ■ • 1 > z k are com “ 
puted Householder transformations, A is the computed Q^Qk-l ■ • • Ql^ z l z 2 ■ • • z k’ 
and B is the computed Q^Qb.i . . . QiBZjZ 2 . . . Z k , the generalized eigenvalue 
problem Ax = \Bx has exactly the same eigenvalues as the problem 

(A + G)x = \(B + H)x where 

l!GI| E ^ |(6.1)k2- t + (45.6)k 2 2" 2t + (92.3)k 3 2- 3t ]llA|| E 

and 

||H|| E s |(6.1)k2- t + (45.6)k 2 2 _2t + (92.3)k 3 2- 3t ]||B|| E 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., October 17, 1974. 
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